X-ray fluorescence holography: going beyond the diffraction limit 
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X-ray fluorescence holography (XFH) is a method for obtaining diffraction-limited images of the 
local atomic structure around a given type of emitter. The reconstructed wave-field represents 
a distorted image of the scatterer electron density distribution; i.e. it is a convolution of the 
charge density distribution with a point spread function characteristic of the measurement. We here 
consider several methods for the iterative deconvolution of such XFH holograms, and via theoretical 
simulations evaluate them from the point of view of going beyond the diffraction limit so as to image 
the electron charge density. Promising results for future applications are found for certain methods, 
and other possible image-enhancement techniques are also discussed. 

PACS numbers: 61.14.-x, 42.40.-i 



I. INTRODUCTION 

X-ray fluorescence holography (XFH) is a promising 
method for directly imaging local atomic structure in an 
element-specific way that haSpbeeji developed experimen- 
tally over the past six yearsoffl. In the first mode of 
measurement, an atom inside an oriented system {e.g. a 
single crystal with low to medium mosaicity) emits a fluo- 
rescent x-ray wave; the interference between the outgoing 
unperturbed reference wave component and the scattered 
object wave component produces a hologramldu. One can 
term this "inside source" holography. A time-reversed 
mode also exists, in which the exciting x-ray wave is the 
reference, and scattered x-rays impinging on a fluorescent 
emitter are the object waves! The emitter is the detec- 
tor of the hologram, and thus this can be called "inner 
detector" holography. The holographic reconstruction in 
either case is achieved by propagating back the wavefield 
from the far field. _ 

Beyond several recent achievements of these methodsQ, 
e.g. in imaging the environment around a dilute semi- 
conductor dopantEl, Jpw-Z atoms in the presence of 
intermediat&j-Z atomsa, and the local environment in a 
quasicrystalQ, another step toward the practical applica- 
tion of holography would be a quantitative link between 
the electronic charge distribution and the holographic 
image. It has already been demonstrated that images 
with accuracies at the diffraction limit can be obtained^. 
But even in this case, the reconstructed wavefield can 
be viewed as a distorted image of the scatterer density 
distribution: i.e. it is the convolution of the charge den- 
sity distribution with a point spread function due to the 
measurement and inversion process. 

When an object is imaged through an optical system, 
information is lost whenever the imaging system cannot 
pass all the spatial frequencies contained in the scene. 
Further loss is produced by aberrations in the optical 
system. The question that we here address is how much 
information can be recovered if we know the point spread 
function (PSF) of the effective optical system associated 



with x-ray fluorescence holography. We will consider here 
several methods for the iterative deconvolution of data so 
as to improve image quality. 

Several procedures for iterative image improvement 
have been developed previously in the fields of astron- 
omy and microscopyfj, and we begin by introducing them 
briefhi Van Cittert developed a very early method in 
192CO; this uses a self-consistent iterative solution to 
the problem. It can be shown that this method is close 
to the gradient search of the maximum likelihood (ML). 
Two other approaches are based on the probability the- 
ory. One developed by LucyEH and RichardsonEj is based 
on Bayes theorem for conditional probability. Another 
one, the maximum entropy method (MEM) seeks the 
'most probable' solution in an under-determined system 
of equationsE3. In this article, we discuss some of these 
methods that seem most appropriate for applying to XFH 
image improvement, quantitatively assess several of them 
via model theoretical calculations, and finally arrive at 
the first deconvolved atomic image at resolution going 
beyond the diffraction limit. 



II. BASIC IMAGING CONCEPTS IN XFH 

We first consider some basic imaging concepts for 
XFH, in particular for an arbitrary atomic charge density 
distribution p (r) . For simplicity, we will also throughout 
this discussion consider only "inner source" holography, 
although the conclusions here can easily be generalized to 
include "inner detector" holography. Neglecting absorp- 
tion or extinction effects (a reasonable approximation if 
we consider only near- neighbor, imaging) , the hologram 
X in the far field k is given bjO: 

X (k) = J T) (k, r) p (r) d 3 r , with 

77(k,r) = 2Re[^e 4 ( fe '- k - r )/(k,r)] ; 



(1) 

here r e is the classical electron radius, r e /(k, r) is the 
scattering factor per electron and the phase factor is 
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given by the path length difference between the wave 
emitted from the origin and the wave scattered by the 
electron located in r. The typical reconstruction algo- 
rithm to calculate the wavefield U in the real space r 

ista 



U(r>) 



{-1} 



(r',k)x(k) d 3 k, with 



,{-1} 



(r',k) = ^Re 



-i(kr'— k-r'1 



(2) 



here a is the measured region in k-space, f2 the volume 
of a, and we have included the scaling factor — to ob- 
tain an estimate of the charge density. Sometimes an 
additional factor is included to compensate the 

'aberration' or distortions produced in the hologram x 
by but we will not explicitly include this here. We 
will for simplicity assume that the scattering factor of an 
electron is isotropic, neglecting the angular dependence 
of the Thompson scattering and near field effects, The 
holographic reconstruction U can thus be considered to 
be a distorted image of the charge density distribution p: 



U(r') 



p(r',r)p(r) d 3 r, with 



M (r',r) = f d 3 k^- 1 >(r',k)r;(k,r); 

J a 



(3) 



with p (r',r) now being the experimental point spread 
function. An additional factor (1 + sign?/ (r',r)) /2 can 
be used to enforce positivity. In summary we have three 
spaces r, k and r' for which we have three functions p (r), 
X (k) and U(r') and the propagators from one space to 
the other are r\ (k, r), r\ (k, r') and p (r',r). 

We now turn to specific methods for iterative decon- 
volution of images. For notational brevity, we in the 
following sections use the symbol • between two func- 
tions to denote integration over the internal variable, e.g. 
f(x, y) ■ g{y) = Jdy f(x, y)g{y). 



III. ITERATIVE IMAGE DECONVOLUTION 

The general problem of the deconvolution is to ob- 
tain the unknown object p (r) from a knowledge of U (r') 
and /i(r', r)_The first method developed by Van Cit- 
tert in 192CO is the most intuitive one and it is based 
on a self-consistent solution. Every step is the difference 
between the reconstructed experimental hologram U (r) 
and the reconstructed simulated hologram [/{»} ( r ) = 
p (r, r')*p{™} (r) generated by a new n t/l -step charge den- 
sity distribution. That is, 

p {n+1} (r) = p<"> (r) + Ap<"> (r) , 
ApM(r) = \u (r) - U {n} (r) 



p {1} (r) = U(r) 



(4) 



Thus, /i(r',r) is in this approximation assumed to be a 
delta function <5(r', r) over each step in the iteration. We 



normally enforce positivity after each step by setting to 
the negative values of the charge density. 

A similar problem is the phase retrieval of a diffraction 
pattern of known intensity from an object with known 
support s (r) . The diffraction pattern provides the am- 
plitude of Fourier transform of the charge density \p\. If 
the object p (r) is constrained within a given support (e.g. 
with the support equal to 1 over the object and outside 
of it), the product of the object p(r) with the support 
function s (r) must be still equal to the object. Therefore 
the complex amplitude of the diffracted wave-field is not 
modified by a convolution with the Fourier transform of 
the support function s: 

FT ~ ~ ~ fK\ 

p = sp > p = s * p . (0) 

If s is sufficiently large ('over-sampling' conditional) 
a certain number of pixels are bounded by this self- 
consistent equation. The main difference with iterative 
deconvolution methods is that in this case, the bigger s 
is, the larger is the number of pixels bounded by (|^), 
and the easier it is to find the solution. If the over- 
sampling condition is satisfied, iterative-approaches such 
as the Gerchberg-Saxtonli3 and FienupEj algorithms can 
be used to solve the phase problem. In the Fourier do- 
main, the Gerchberg-Saxton algorithm can be expressed 
as: 



p{n+l} 



Pmeas e 



phase{p ( " } *s} 



where the n and (n + 1) estimates for density are 
given, and |/5meas| is the measured amplitude of the 
Fourier transform of the object. 

Another iterativ e . ap proach is given by the Lucy- 
Richardson methodtlHU. This method is based on a self- 
consistent iterative solution of the Bayes theorem for con- 
ditional probabilities, considering p (r', r) as a probabil- 
ity that r' will fall in the interval (r', r' + dr') when it is 
known that r = r'. A self-consistent iterative approach 
can in this case be expressed astil: 



p {n+l} ( r ) = 



p {n} ( r ) 
V 1 ' 1 I /j(r',r")*p(">(r") 



(6) 



with p T (r,r') being the transpose of p. 



IV. MAXIMUM LIKELIHOOD AND ENTROPY 

In another set of methods, we want to minimize the dif- 
ference between the measured hologram in the far field 
X (k) and the hologram generated by a given charge den- 
sity distribution p^ (r): 



Q 



d 3 k 



X(k)--)(k,r)-p w (r) 



(7) 



where the unknowns are p^ (rj) , i — 1...N over some 
grid spanning the object. The gradient can be computed 
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VQ = 2rf (r, k) • 7] (k, r 



■ P {n} (r) 



X(r) 



(8) 



with if being the transpose conjugate of rj. The Hessian 
of Q is in this case nothing else but the point spread 
function and is given by: 



H 



d 2 Q 



dpi' 1 } (n) V™ } (r 2 ) 
2rf (ri,k) • i] (k,r 2 ) . 



(9) 



The single pixel appr oximat iorr tt H consists in assuming 
that H is diagonal, so that the steps will be: 



Ap { " } (r) = -J-VQ r 

tlr.r 



(10) 



Under this approximation, the gradient search is equiva- 
lent to the Van Cittert method, since H T r ~ (^) (see 
proof in Appendix |A[) and rjl -1 } = However, for 

single energy holograms, at the inversion-symmetric po- 
sition to the real atoms, there will appear a twin or ghost 
image. Therefore one should consider at least the sym- 
metric term H r _ r : 



H r , r ' 

Apf™> (r) 



H r , r ' (Sr',r + #r',-r) 



\ VQ r 



(11) 

(.12) 



An additional consideration here however, is that at ev- 
ery energy, for a given position, and in case of a point 
scattering, the twin image cancels the real image, so that 
a divergence arises in Eq. ([l^). This can be avoided 
by introducing an anisotropy in the scattering factor h£=. 
tween the forward and backward directions. It is knownEHI 
that, if we use a small Gaussian charge distribution for 
every point instead of point scatterers, the scattering 
becomes stronger in the forward direction, and iterative 
algorithms converge more quickly. Conjugate gradient 
methods do not require the calculation of the Hessian, 
but partial knowledge of it is helpful. Extra information 
can also be used to speed up the process, and for ex- 
ample, we know that the scattering charge distribution 
is real and positive. This can be enforced by taking the 
real part at every step and by setting to p n = when 
p n (r) < 0. A further piece of information is the total 
number of scattering electrons Fo, and we can use La- 
grange multipliers to perform a constrained best fit via: 



Q + p (F - F ) , with 
F = 5>( ri ). 



(13) 



Here the minimization is performed adjusting the multi- 
plier at every step: 



A <°r = -H^,(VQ r ,+/3VF r 0, 

[F - F + (VF • fi- 1 • VQ)] 



A 



VF • H- 1 • VF 



(14) 



The Lagrange multiplier method keeps the number of 
electron fixed, however after we apply each step, we en- 
force positivity in the charge density by setting to the 
negative values of the charge density, therefore the actual 
number of electrons will tend to grow at every iteration. 
In the actual implementation of the algorithm, we arti- 
ficially reduce the step size by a factor of 1/5 to avoid 
cutting too much in the negative density regions. 

The deconvolution of spectra with an incomplete 
Fourier image is_also often treated with the maximum 
entropy methodtJ. When the number of unknowns p (r 4 ) 
is superior to the number of equations, it is necessary to 
include further a-posteriori information. As more than 
one solution can satisfy the set of equations, we seek the 
most probable solution, i.e. the one with the maximum 
entropy. In practice we want to minimize the functional: 



S +A'Q, 

S =y\p(r i )ln / o(ri) 

' " 1. 



(15) 



where A' is a Lagrange multiplier. The procedure is sim- 
ilar to the previous one, with two main differences: (i) 
MEM enforces positivity, so that a pure gradient min- 
imization could lead to unphysical values that need to 
be chopped, (ii) S is highly nonlinear, and simple gradi- 
ent search is not efficient. An interesting property of S 
is that its Hessian is diagonal; therefore the single pixel 
approximation can be usedEll. 



V. NUMERICAL SIMULATIONS 

We have initially tried all the above mentioned decon- 
volution methods on a test system consisting of a charge 
distribution described by a 3D matrix of 50 3 points with 
0.2 A spacing. This matrix was set to except for two 
Gaussian charge distributions centered in the x-z plane 
representing two atoms separated by 3.5 A as shown in 
Fig. (0a). Figures (0b-d) now show reconstructed images 
via several different methods obtained from a hologram 
(Fig. |^a) defined at an energy of 17 keV with the points 
distributed over a solid angle with azimuth varying from 
0° to 360° and polar angle from 0° (parallel to the z- 
axis) to 80°. The reconstruction in Fig. (|]f) has been 
obtained using both holograms at 17 keV and 18 keV 
(Fig. |]b). The standard direct Barton method in Fig. 
(0b) shows a weak twin image, and some distortion of the 
imaged charge distributions, with the difference between 
Fig. (0a) and Fig. (0b) essentially being the PSF. The 
number of electrons estimated from the Barton method 
is one order of magnitude bigger than the original one. 
This could be ascribed to the fact that many of these 
electrons are not contributing to the hologram, as they 
cancel each other. In trials with the Van Cittert and ML 
(Fig. 0c) methods using the double pixel approximation, 
we have found that the convergence is very good: less 
then 10 iterations are necessary for convergence. This is 
due to the fact that the holographic reconstruction is a 
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FIG. 1: real space image of (a) the original charge density dis- 
tribution - Fo = 50, (b) Barton reconstruction - F = 400, (c) 
Maximum Likelihood - F = 550, (d) Maximum Entropy and 

- F = 300, (e) constrained ML - F = 60, (f) Constrained ML 

- F = 60. (b-e) were obtained from the hologram at 17 keV 
in Fig. (^a) while (f) was obtained from the simulated holo- 
grams at 17 and 18 keV shown in Fig. (^,b). The colormaps 
are scaled for each picture. 



relatively good image of the charge density. However the 
number of electrons used to fit the data was one order of 
magnitude bigger than that of the original image. An- 
other relevant feature of this method is the fact that the 
image becomes more structured. This is due to the di- 
vergence in the inverse of the Hessian using single energy 
holograms when H r r = H r r . 

Tests using the Lucy-Richardson method have shown 
bad convergence; further studies of this method need to 
be done, but one of the reasons could be some divergence 
caused by the denominator in equation (JsJ) . 

At this point we turned to the maximum entropy 
method. As we had already developed the ML algo- 
rithm with single and double pixel approximations, we 
could easily apply the algorithm developed by Cornwell 



FIG. 2: Top-Simulated holograms from two Gaussian charge 
distributions as shown in Fig. (0a) at 17 keV and (jjjb) 18 keV 
after smoothing. Center- Holograms generated by the fitted 
charge density distribution at the same energies (lllc) 17 keV 
and (firl) 18 keV. Bottom- Difference between top and center: 
[J = - |c and |l|f = §b - |d. 



and EvansciJ so as to employ the maximum entropy cri- 
terion. This algorithm tries to maximize the entropy S 
while minimizing the likelihood Q with a fixed number 
of electrons F using Lagrange multipliers. Several tests 
have shown that this algorithm (at least the one imple- 
mented by us) was not able to minimize Q, S and F at 
the same time; the result is shown in Fig. (|lji). How- 
ever, we did find that the maximum likelihood and the 
minimum number of electrons produces a good quality 
reconstruction, with somewhat improved resolution, and 
a reasonably good FtJ; these results are shown in Fig. 
(|]e) and they can be compared with those in (|l]b). The 
number of electrons is only 20% bigger than the original 
one. 

The single energy iterative reconstruction shows how 
the divergence in the Hessian due to the twin image prob- 
lem produces several artifacts, which are removed when 
using slightly different energies (Fig. |l]f). The reason 
why the number of electrons does not correspond to the 
original Fo even when applying the constraint is that af- 
ter the steps are applied, we impose positivity by setting 
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FIG. 3: Top- Original charge distribution (F = 1.8 x 10 4 ) - 
(a) XZ plane, and (b) 3D view. Center- Standard holographic 
reconstruction (F = 7 x 10 4 ), (c) and (d) the same views as 
(a) and (b), respectively. Bottom-Deconvolved charge density 
distribution (F = 2.9 x 10 4 ), (e) and (f) show the same 
views as (a) and (b). The colorbar shows the maximum and 
minimum value. The 3D images are obtained by maximum 
voxel (volume-pixels) projection: the 3D matrix is rotated 
with respect to a plane (screen). On every pixel of the screen 
we project the maximum value of the voxels on top of the 
pixel itself. 
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FIG. 4: Top-Simulated holograms from a cluster of 1000 
atoms at (a) 8 keV and (b) 9 keV after smoothing. Center- 
Holograms generated by the fitted charge density distribu- 
tion at the same energies (a) 8 keV and (b) 9 keV. Bottom- 
Difference between top and center: e = a - c and f = b - d 



to the negative values of the density. 

Finally, we applied the constrained ML method to the 
same case, using a double-energy data set. This result 
is shown in Fig. (|l]f) and represents the most significant 
enhancement in image quality relative to the standard 
method. 

To explore the constrained ML algorithm further, we 
carried out simulations based on a much larger and more 
realistic cluster of atoms representing CdTe; the charge 
density distribution in the neighborhood of the origin 
is shown in Fig. ([|a,b). We tried to keep a minimum 
amount of information to show the robustness of the 
method: We simulated the hologram produced by a CdTe 
cluster of 1000 atoms at only two energies of 8 and 9 keV. 
At these two energies the holographic reconstruction us- 
ing the Barton method shows a poor image quality as 
shown in Fig. (||c,d). We also introduced a randomly 
distributed noise level of 10~ 3 on every data point (cor- 



responding to 10 6 counts in each point), with the points 
distributed over a solid angle with azimuth varying from 
to 360 degrees and polar angle from (perpendicular 
to the [001] surface normal direction) to 80 degrees, with 
a 1-degree step size in each direction. Since we are inter- 
ested in imaging the first neighbors, we finally carried out 
all imaging with a smoothed hologram based on 2 degree 
smoothing in every direction. Although the cluster size 
is several nanometers, we will only try to reconstruct up 
to 8 A setting to anything outside this range. 

The holograms generated by this charge after smooth- 
ing are shown in Fig. (0a) at 8 keV and Fig. (J4|d) at 
9 keV. We again choose as a first reference step and set 
of images the standard holographic reconstruction in Fig. 

,b). Comparison between the standard reconstruction 
and the deconvolved images are shown in Fig. (D^f). 

The number of electrons within the reconstructed re- 
gion is F = 1.8 x 10 4 in the original picture, F = 7 x 10 4 
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FIG. 5: behavior of the fit quality (Q) normalized to the vol- 
ume in k-space, the number of electrons (F) and the Lagrange 
multiplier A. The fit is obtained quickly, while F decreases 
much more slowly. The multiplier follows the behavior of F 
as it tries to bring it back to Fq ~ 10 4 . 



in the standard reconstruction, and F = 2.9 x 10 4 in the 
deconvolved image. Thus, the deconvolved image is a 
factor of 1.5 larger than the original image in number 
of electrons, but a factor of 2 better than the standard 
reconstruction. This result is worse than what we ob- 
served for the two Gaussian charges, and can be ascribed 
to the difficulty to fit point like structures (atoms) at 
the diffraction limit scale instead of smooth distributions. 
The peak heights of the charge density in the three im- 
ages varies somewhat differently, the maximum number 
of electrons in one voxel in the original image is 50, 0.38 
in the standard reconstruction, and 14 in the deconvolved 
image, as shown more quantitatively in the colorbars of 
Fig. (|J). The reason for this is that in the original im- 
age, every electron of one atom is concentrated in a single 
voxel, while the standard reconstruction has many elec- 
trons spread over many more voxels. The deconvolved 
image has a total of fewer electrons than the standard 
image, but much higher peaks, due to the enhanced res- 
olution relative to the standard reconstruction. Com- 
paring image quality between the standard (Fig. ||c,d) 
and deconvolved (Fig. ||e,f) approaches, we see a definite 
enhancement of resolution in the deconvolved images, 
even though some artifacts remain that are of compa- 
rable strength to the actual charge-density images. The 
holograms generated by the cluster in Fig. (^a,b) are 
well described in the low frequency range by the fitted 
holograms in Fig. (|c) (8 keV) and Fig. (|d). bmce m 
the iterations we only fitted the charge density within the 
first 8 A setting the distribution to outside this range, 
only the low frequency components where fitted, as the 
residual difference shows in Fig. (^e-f). 

In Fig. (||), we present some results dealing with the 
rapidity and type of convergence observed. The maxi- 



mum likelihood is obtained within the first 10-20 iter- 
ations (Fig. |^-top). The total number of electrons is 
slowly reduced in the following ~ 500 iterations (Fig. ||- 
center). The Lagrange multiplier A follows the behavior 
of the constraint: as F changes from the desired level Fo, 
A varies in order to bring F back to the desired level. 



VI. CONCLUSIONS AND FUTURE 
POSSIBILITIES 

In conclusion, we have demonstrated that iterative de- 
convolution can be applied to XFH imaging, thus im- 
proving the quality of the images beyond the diffraction 
limit, and that several levels of approximation can be 
used. From our test cases, the constrained maximum 
likelihood method appears to be the most promising. In 
these first theoretical simulations related to image im- 
provement, we have approximated the scattering factor 
as isotropic, thus neglecting the angular dependence of 
Thomson scattering and near-field effects. As a compu- 
tationally relevant issue, we have also approximated the 
Hessian to be a sparse matrix where only the two diag- 
onals are non-zero. All of these approximations could 
be relaxed in future applications, including terms close 
to the diagonals, with of course additional complexity 
then being incurred in the deconvolution algorithm and 
the length of time it takes. For our calculations, the 
storage and calculation of the full Hessian of 80 6 points 
required considerable memory, but computers capable of 
such calculations should be available as desktops in a few 
years. The gradient was also calculated neglecting the 
angular dependence of the scattering factor. Using these 
approximations, every iteration involved the calculation 
of a hologram from a matrix of 80 3 voxels, the recon- 
struction and required about 5 minutes on a Pentium II 
computer, for a total of one night of calculation. Our 
image deconvolution has been obtained by enforcing a 
real and positive charge density distribution and limited 
number of electrons. However, looking ahead to further 
image improvement, once the image is deconvolved, one 
could make use of 'atomicity' in the next step by placing 
atoms in the location of the peaks, and then perform- 
ing an optimization using the full angular dependence 
of the scattering factor, while simultaneously adjusting 
the atomic number and position of these atoms. This 
study thus suggests several fruitful directions for further 
exploration and exploitation of image deconvolution in 
x-ray fluorescence holography that should considerably 
enhance its power for structural studies. 
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APPENDIX A: APPENDIX-CALCULATION OF 
THE HESSIAN AND THE GRADIENT 

The Hessian does not need to be known in great detail. 
It can be pre-calculated before starting the iterative al- 
gorithm neglecting the angular dependence of Thompson 

scattering factor ^ = 1+cos 2 ek,r and using the follow- 
ing approximation f cos 2 (kr — k • r) k 2 d 2 k ~ ifij. and 

a(k) 

J sin 2 (k • r) k 2 d 2 k ~ §n fe , with Q k = / , fc . k 2 d 2 k, k is 

a(k) 

a unit vector, and u(k) is the region where the hologram 
is measured at a given energy. This yields: 

H r r = V Ak [ 4 cos 2 (kr - k ■ r) k 2 d 2 k, 

k J<j(k) r " 



H r r = V Ak [ -^cos(kr-k 



-kr) 



x cos (kr + k • r) k 2 d 2 k 

% Ek Afc/ a(fc) [cos 2 (kr) - sin 2 (k ■ r)] fc 2 d 2 k, 
£ Ek [cos 2 (kr) - i] n k Ak 



2 ^k 



Y^ k cos (2kr) O fc Ak ; 



(Al) 



The gradient is calculated by simulating the hologram 
t^W (k) generated by a charge density distribution 
p{n} ( r ) 5 anc j propagating back the hologram to the real 

space. Using U (r) = (jr^j U (r), eq. (||) can be written 



as: 



VQ = 2 [u {n} (r) - U (r) 



[/{«} 



x {n} (k) 



7?t (r,k)- X { " } (k) 

-^Re / e -*(*»--kT) {«} (k)d 3 k; 

r?(k,r).pW (r) 

2R c y ^ e i ( fcr - k - r V {n} 0)d 3 r . (A2) 



These calculations can be very long as we need to sum 
over r(x,y,z) for every k(k,$, (f). Supposing that the 
number of points in x,y, z,ft,(p is N ~ 100 for every 
variable, we need to calculate the product in the integral 
of eq. (A2) for every possible value, i.e. N 5 products and 
sums. The computation of x^ s ) can be performed 
faster if we separate the sum over the variables x, y, z. 

" (r) Ar 3 ; 



We begin by calculating ^ n > (r) 



{"} ( r \ — JLf,ikr n {n} 



we can express Eq. (A2) as: 



Xn (k) ~2R£' 



-ik-r 



C {n} (r) 



(A3) 



which is a simple Fourier transform. However the holo- 
grams are usually measured in spherical coordinates, and 
to obtain the values of Xn (k) in the measured positions, 
we either need to perform an interpolation, or a convo- 
lution. Another approach is to use directly the spherical 
coordinates in the kernel: 



gik-r gikx sin $ cos <p ^iky sin •& sin ip^ikz cos $ 



(A4) 



and Eq. (A3) becomes: 

X {n} (k,tf,v) 2ReV Ak.x.v^Y] B k , y ,i>, v 

xV C k ^^(x,y,z). (A5) 

* — ' z 

Each sum requires N A calculation, reducing the time to 
perform the calculation by ~N/3. The same trick can be 
applied for the holographic reconstruction. 
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